Depletion potentials in highly size-asymmetric binary hard-sphere mixtures: 
Comparison of accurate simulation results with theory 
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We report a detailed study, using state-of-the-art simulation and theoretical methods, of the ef- 
fective (depletion) potential between a pair of big hard spheres immersed in a reservoir of much 
smaller hard spheres, the size disparity being measured by the ratio of diameters q = cr s /<?(,. Small 
particles are treated grand canonically, their influence being parameterized in terms of their packing 
fraction in the reservoir, rf s . Two specialized Monte Carlo simulation schemes -the geometrical 
cluster algorithm, and staged particle insertion- are deployed to obtain accurate depletion poten- 
tials for a number of combinations of q < 0.1 and if s . After applying corrections for simulation 
finite-size effects, the depletion potentials are compared with the prediction of new density func- 
tional theory (DFT) calculations based on the insertion trick using the Rosenfeld functional and 
several subsequent modifications. While agreement between the DFT and simulation is generally 
good, significant discrepancies are evident at the largest reservoir packing fraction accessible to 
our simulation methods, namely rf s — 0.35. These discrepancies are, however, small compared to 
those between simulation and the much poorer predictions of the Derjaguin approximation at this 
rf s . The recently proposed morphometric approximation performs better than Derjaguin but is 
somewhat poorer than DFT for the size ratios and small sphere packing fractions that we consider. 
The effective potentials from simulation, DFT and the morphometric approximation were used to 
compute the second virial coefficient B2 as a function of rf a . Comparison of the results enables an 
assessment of the extent to which DFT can be expected to correctly predict the propensity towards 
fluid fluid phase separation in additive binary hard sphere mixtures with q < 0.1. In all, the new 
simulation results provide the first fully quantitative benchmark for assessing the relative accuracy 
of theoretical approaches for calculating depletion potentials in highly size-asymmetric mixtures. 



I. INTRODUCTION 

Much of condensed matter physics and chemistry is 
concerned with simplifying the description of a complex 
many body system by integrating out certain subsets 
of the degrees of freedom of the full system. Thus in 
treating atomic and molecular solids and liquids one of- 
ten resorts to integrating out, usually approximately, the 
higher energy quantum mechanical degrees of freedom of 
the electrons in order to obtain an effective interatomic 
or intermolecular potential energy function that is then 
employed in a classical statistical mechanical treatment 
to study the properties of condensed phases. Similarly in 
metallic systems integrating out the degrees of freedom 
of the conduction electrons leads to an effective Hamil- 
tonian for the screened ions or pseudo-atoms. When one 
turns to complex, multi-component fluids such as col- 
loidal suspensions or polymers in solution the basic idea 
is similar: one integrates out the degrees of freedom of 
the small species to obtain an effective Hamiltonian for 
the biggest species; for an illuminating review see [1]. In 
this case all species can be treated classically and the 
formalism is essentially the famous one of McMillan and 
Mayer [2] who developed a general theory for the equi- 
librium properties of solutions. These authors and many 
subsequently recognized that integrating out is best per- 
formed when the smaller species (those constituting the 
solvent) are treated grand canonically. Obtaining the full 
effective Hamiltonian is a tall order. A first step in any 



theoretical treatment is to determine the effective poten- 
tial between a single pair of the biggest particles in a sea 
of the smaller species. There is a long history of work 
in this field. Different statistical mechanical techniques 
have been developed to calculate these potentials for dif- 
ferent types of complex fluid. Well-known examples of ef- 
fective two-body potentials, forming cornerstones of col- 
loid science, are the DLVO potential for charged colloids 
and the Asakura-Oosawa depletion potential for colloid- 
polymer mixtures. Further examples, including polymer 
systems, are given in [1, 3]. 

The problem of determining the effective two-body in- 
teraction is particularly challenging to theory and com- 
puter simulation in the situation where - and we special- 
ize to a binary fluid - the bigger particles are much larger 
than the smaller ones. Sophisticated theoretical and sim- 
ulation techniques, that will provide accurate results for 
the effective potential, are only now becoming available. 

Our present focus is on a simple model for a suspension 
of a binary mixture of big and small colloidal particles, 
both species suitably sterically stabilized, i.e. we consider 
a highly size asymmetric binary mixture of hard spheres. 
This can also serve as a crude model of a mixture of 
colloids and non-adsorbing polymer and can be regarded 
as a reference system for a mixture of size asymmetric 
simple fluids. The hard sphere mixture is important in 
that it provides a testing ground for theories of the effec- 
tive potential between two big hard spheres immersed in 
reservoirs of small ones, with different packing fractions 
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rf sl and indeed for theories of the full effective Hamilto- 
nian obtained when the degrees of freedom of the small 
species are integrated out fully [4]. Generally speaking, 
the more asymmetric the mixture the more difficult it is 
to treat both species on equal footing and the more nec- 
essary it is to perform some integrating out to obtain an 
(approximate) effective Hamiltonian. Often this proce- 
dure is only feasible at the level of a Hamiltonian consist- 
ing of a sum of (effective) pair potentials, as obtained by 
considering a single pair of the big particles, along with 
zero and one body terms [1, 4]. One of the advantages 
of the hard sphere system is that geometrical considera- 
tions indicate that three, four etc body terms become less 
important as the size ratio q = o s j(Jb becomes small. a s 
and <7b denote the diameters of the small and big species, 
respectively. Thus provided one can calculate an accu- 
rate effective pair potential, the pair description alone 
determines an effective Hamiltonian that should provide 
an excellent description of the big-big correlation func- 
tions and the phase behaviour of the binary hard sphere 
mixture when q is sufficiently small [4] . Note that in this 
paper we consider additive hard sphere mixtures so that 
the big-small diameter cr^s = (<T(, +<t s )/2. 

Since the studies of the phase behaviour of the hard 
sphere mixture by Dijkstra et al [4], whose simulations of 
an effective one-component system used a rather crude 
approximate pair potential, there have been several new 
developments in the theory of effective potentials. Most 
of these are based on Density Functional Theory (DFT). 
It is important to assess whether i) the potentials derived 
in recent studies are accurate and ii) use of these might 
lead to different predictions for the properties of the mix- 
ture. In order to make such assessments it is necessary 
to have accurate simulation results for the effective po- 
tential. Employing state of the art techniques we provide 
what we believe are the most accurate results currently 
available for q < 0.1 and packing fractions rf s up to 0.35 
and make direct comparisons with the results of theoreti- 
cal approaches. The simulation techniques we employ do 
not allow us to work at very high values of rf s but they 
do allow us to compute accurate effective potentials, for 
a range of size ratios, in the regime of small sphere pack- 
ing fractions where the putative (metastable) fluid-fluid 
phase separation is predicted to occur [4]. By calculat- 
ing the second virial coefficient associated with the effec- 
tive potential we make new estimates of the value of r] r s 
where the onset of this elusive phase transition occurs. Of 
course real colloidal systems may not reach equilibrium 
on experimental timescales, particularly if the effective 
potential exhibits significant repulsive barriers [5]. Nev- 
ertheless, knowledge of the underlying phase behaviour is 
important for interpreting dynamical observations, such 
as whether gelation or glassiness might set in [6, 7]. 

The comparison between DFT results and simulation 
addresses recent suggestions [8, 9] that no existing theo- 
retical framework is reliable for calculating effective po- 
tentials for small values of q and physically relevant val- 
ues of rf s . We examine and refute these suggestions in 



the light of our present results. 

The effective potential between two big hard spheres 
takes the form: 

4>eff{rbb) = 4>bb{nb) + W(r bb ) (1) 

where 4>bb is the bare hard sphere potential between two 
big spheres and W is the so-called depletion potential. 
This is attractive for small separations r b b of the big 
spheres but decays in an exponentially damped oscilla- 
tory fashion at large separations. The physics of the at- 
traction is well understood: the exclusion or depletion of 
the small spheres as the big ones come close together re- 
sults in an increase in free volume available to the small 
species leading to an increase of entropy. If this attrac- 
tion is sufficiently strong it can give rise to fluid-fluid 
phase separation. Such a phase transition is driven by 
purely entropic effects: recall that all the bare interac- 
tions in the mixture are those of hard spheres. Of course 
the concept of an attractive depletion potential between 
colloids dispersed in a solution of non-absorbing polymer, 
or other depletants has a long history. The recent book 
by Lekkerkerker and Tuinier [10] describes this and the 
general importance of depletion interactions in colloidal 
systems. 

We have emphasized that the effective pair potential 
is a key ingredient in an effective Hamiltonian descrip- 
tion of the mixture. However, this object is also impor- 
tant in its own right since it can be measured experi- 
mentally for colloidal systems using various techniques. 
More specifically, the effective potential between a single 
colloid, immersed in a suspension of small colloidal par- 
ticles or non-adsorbing polymer, and a flat substrate has 
been measured; see for example [11] and the comparisons 
made between DFT results and experiment [12]. Crocker 
et al [13] measured the effective potential between two big 
PMMA particles immersed in a sea of small polystyrene 
spheres and observed damped oscillations at high small 
sphere packing fractions. Subsequently comparisons were 
made with DFT results [14]. Ref. [15] provides a recent 
review of direct experimental measurements of effective 
interactions in colloid-polymer systems. 

A. Previous Simulation Studies 

In general, the task of accurately measuring effective 
potentials in highly size-asymmetrical fluid mixtures us- 
ing traditional simulation methods such as molecular dy- 
namics (MD) or basic Monte Carlo (MC) is an extremely 
challenging one. The difficulty stems from the slow re- 
laxation of the big particles caused by the presence of 
the small ones. Specifically, in order for a big particle 
to relax, it must move a distance of order its own diam- 
eter a b - However, for small size ratios q, and even at 
quite low values of if s , very many small particles are typ- 
ically to be found surrounding a big particle and these 
hem it in, greatly hindering its movement. In computa- 
tional terms this mandates a very small MD timestep in 
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order to control integration errors, while in MC, a very 
small trial step-size must be used in order to maintain 
a reasonable acceptance rate. Consequently, the com- 
putational investment required to simulate highly size 
asymmetric mixtures by traditional means is (generally 
speaking) prohibitive at all but the lowest packing frac- 
tions of small particles. 

Owing to these difficulties, most previous simulation 
studies of hard sphere mixtures [8, 9, 16-18] have adopted 
an indirect route to measuring effective potentials in the 
highly size asymmetrical regime q < 0.1 based on mea- 
surements of interparticle force. The strategy rests on 
the observation that the force between two big parti- 
cles can be expressed in terms of the contact density of 
small particles at the surface of the big ones [17, 19]. By 
measuring this (angularly dependent) contact density for 
fixed separation r bb of the big particles and repeating for 
separations ranging from the minimum value r b b = <?b to 
rbb = oo, one obtains the force profile F(rb b ). This can in 
turn be integrated to yield an estimate of the depletion 
potential. Generally speaking, however, the statistical 
quality of the data obtained via this route seems typi- 
cally quite low, particularly at small q and high nf s . This 
presumably reflects the difficulties of measuring contact 
densities accurately and the errors inherent in numerical 
integration. 

Only a few studies have tried to measure the deple- 
tion potential directly for q < 0.1 (see ref. [20] for a hard 
sphere study and refs. [21, 22] for more general poten- 
tials). In common with the present work, these studies 
deployed a cluster algorithm (to be described in sec. II A) 
to deal efficiently with the problem of slow relaxation out- 
lined above. However, they treated the small particles 
canonically rather than grand canonically, which compli- 
cates comparison with theoretical predictions which are 
typically formulated in terms of an infinite reservoir of 
small particles [23]. Furthermore it seems that no pre- 
vious simulation studies have discussed (in any detail) 
finite-size effects in measurements of depiction potentials, 
the role of which we believe to be particularly significant 
at large size asymmetries. Consequently, while previous 
work has evidenced good qualitative agreement between 
simulation and theory, there is to date a lack of data from 
which one can make confident comparisons between the 
various theoretical approaches. This is remedied in the 
present work. 



B. Previous Theoretical Studies 

There are many of these and they are based on a vari- 
ety of techniques. Integral equation treatments abound 
and these are summarized nicely in the recent article by 
Bol^an et al [24] . A related rational function approach was 
used recently by Yuste et al [25]. Tackling highly asym- 
metric mixtures via integral equation methods, where one 
treats both species on equal footing, is notoriously diffi- 
cult and making systematic assessments of the reliability 



of closure approximations is not straightforward and, of 
course, requires accurate simulation data. 

Density functional (DFT) treatments are arguably 
much more powerful. There are several different ways 
of calculating the effective (depletion) potential between 
two big hard spheres in a reservoir of small hard spheres 
or more generally of calculating the effective potential be- 
tween two big particles in a reservoir of small ones with 
arbitrary interactions between bb, bs and ss. The first 
method is to fix the centres of the big (b) particles a dis- 
tance rb b apart and then compute the grand potential 
of the small (s) particles in the external field of the two 
fixed b particles as a function of r b b for a given size ratio 
q and a reservoir packing fraction of rf s . This method re- 
quires only a DFT for a single component fluid, the small 
particles. The big particles are fixed so they simply ex- 
ert an external potential on the small ones. DFT for 
one-component hard spheres is very well-developed; very 
accurate functionals exist and these are suitable for treat- 
ing the extreme inhomogeneities that arise for small size 
ratios q. It follows that this brute force method should be 
rather accurate. Its drawback is that the density profile 
of the small particles has cylindrical symmetry requir- 
ing numerically accurate minimization of the free energy 
functional on a two-dimensional grid. 

Goulding [26, 27] performed pioneering brute force cal- 
culations using the Rosenfeld fundamental measure the- 
ory (FMT) [28] for a system with q = 0.2 and pack- 
ing fractions rf s up to 0.314. This method has been 
refined recently by Bo^an et al [24] who employed vari- 
ous hard-sphere functionals for more asymmetric systems 
and higher values of -q r s . These authors, see also Oettel 
et al [29], also used DFT to calculate the depletion force 
directly using the formula due to Attard [17, 19] that re- 
lates the force to the density profile of small spheres in 
contact with a big sphere. Once again the density profile 
has cylindrical symmetry and careful numerical methods 
are required. 

A popular DFT method for hard-sphere systems is 
based on what has become known as the insertion trick 
or insertion method [14, 18]. This is a general procedure, 
see Sec. Ill A, for calculating the depletion potential be- 
tween a big particle and a fixed object, e.g. a wall or 
another big particle, immersed in a sea of small particles. 
The advantage of the method is that one requires only 
the equilibrium density profile of the small species in the 
external field of the isolated fixed object and this profile 
clearly has the symmetry of the single fixed object. For 
the case of two big spheres the profile p s (r) has spherical 
symmetry. The disadvantage is that the theory requires 
a DFT for an asymmetric mixture, albeit in the limit 
where the density pb of the big particles approaches zero: 
Pb — > 0. For hard-sphere mixtures the insertion method 
is straightforward to implement and [14] provides a series 
of comparisons, using the Rosenfeld FMT [28], with the 
simulation data that existed in 2000. 

Further comparisons between results of the DFT in- 
sertion method and simulation studies were made in 
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Refs. [8, 9, 24, 29]. In the present paper we seek to make 
more quantitative comparisons, taking into consideration 
the improved accuracy of our new simulation results and 
the availability of improved DFTs for hard sphere fluids. 

There is a further theoretical approach to calculating 
depletion potentials developed very recently in Ref. [29] 
and employed subsequently in Ref. [24]. This approach 
is based on morphological (or morphometric) thermody- 
namics [30] . The basic idea is that the depletion potential 
is (essentially) the solvation free energy of the dumbbell 
formed by the two big spheres and that this quantity can 
be separated into geometrical measures, namely the vol- 
ume, surface area and integrated mean and Gaussian cur- 
vatures. The coefficients of these measures are geometry 
independent thermodynamic quantities, i.e. the pressure, 
the planar surface tension and two bending rigidities all 
of which can be obtained from simulations or from DFT 
calculations of the single component fluid performed for 
a simple geometry. 

The paper is arranged as follows: Section II describes 
the grand canonical simulation methods that we have em- 
ployed for determining the depletion potential between 
two big spheres for size ratios q from 0.1 to 0.01. In Sec- 
tion III we summarize briefly the DFT insertion method 
and the three hard sphere functionals that we employ in 
our present calculations. We also discuss some of the lim- 
itations of the parameterization of the depletion potential 
introduced in [14]. Some details of the morphometric ap- 
proaches are also given here. Results are presented in 
Section IV. As a test of our simulation method we de- 
termine the depletion potential for two big hard spheres 
in a solvent of non-interacting point particles that have 
a hard interaction with the big spheres. For this case the 
depletion potential is known analytically-it is the ven- 
erable Asakura-Oosawa potential [31, 32] of colloid sci- 
ence. For the additive binary hard sphere case we make 
comparisons between results of simulation, DFT inser- 
tion method, the morphometric approach and the Der- 
jaguin approximation [33] for the depletion potential. We 
also compare simulation, DFT and morphometric results 
for the associated second virial coefficient B 2 (r]g). The 
latter provides a valuable indicator of the propensity of 
the bulk binary mixture to phase separate into two fluid 
phases [34, 35]. We conclude in Section V with a discus- 
sion. 



II. SIMULATION METHODS 

A. Geometrical cluster algorithm 

An efficient cluster algorithm capable of dealing with 
hard spheres mixtures was introduced by Dress and 
Krauth in 1995 [36]. It was subsequently generalized 
to arbitrary interaction potentials by Liu and Luijten 
[21, 37] who dubbed their method the Geometrical Clus- 
ter Algorithm (GCA). Here we describe the application 
of the GCA to a size asymmetrical binary mixture of hard 



spheres. 

The particles comprising the system are assumed to 
be contained in a periodically replicated cubic simula- 
tion box of volume V. The configuration space of these 
particles is explored via cluster updates, in which a sub- 
set of the particles known as the "cluster" is displaced via 
a point reflection operation in a randomly chosen pivot 
point. The cluster generally comprises both big and small 
particles and by virtue of the symmetry of the point re- 
flection, members of the cluster retain their relative posi- 
tions under the cluster move. Importantly, cluster moves 
are rejection-free even for arbitrary interparticle interac- 
tions [21]. This is because the manner in which a cluster 
is built ensures that the new configuration is automati- 
cally Boltzmann distributed. 

For hard spheres, the cluster is constructed as follows: 
one of the particles is chosen at random to be the seed 
particle of the cluster. This particle is point-reflected 
with respect to the pivot from its original position to 
a new position. However, in its new position, the seed 
particle may overlap with other particles. The identi- 
ties of all such overlapping particles are recorded in a 
list or "stack" . One then takes the top-most particle off 
the stack, and reflects its position with respect to the 
pivot. Any particles which overlap with this particle at 
its destination site are then added to the bottom of the 
stack. This process is repeated iteratively until the stack 
is empty and there are no more overlaps. 

In this work we shall be concerned with measurements 
of the radial distribution function gbb(r) for a system con- 
taining a pair of big particles among many small ones. To 
effect this measurement we modify the GCA slightly as 
follows: we choose one big particle to be the seed particle, 
which we place randomly within a shell abb < r < L/2, 
centred on the second big particle, with L the linear box 
dimension. The location of the pivot is then inferred from 
the old and new positions of the seed particle. Thereafter 
clusters are built in the standard way. This strategy en- 
sures that we efficiently sample separations of the big 
particles that lie in the range for which <7f>&(r) can sensi- 
bly be defined for hard spheres in a cubic box. 

Small particles are treated grand canonically in our 
simulations. In practical terms this means that in paral- 
lel with the cluster moves, we implement insertions and 
deletions of small particles, subject to a Metropolis ac- 
ceptance criterion governed by an imposed chemical po- 
tential. The choice of chemical potential controls the 
reservoir packing fraction of small particles. 

For the systems of interest in this work, we find that 
the GCA is efficient for reservoir packing fractions rf s < 
0.2. Above this value one finds that practically all the 
particles join the cluster, which merely results in a trivial 
point reflection of the entire system. For single compo- 
nent fluids this problem can be ameliorated by biasing 
the choice of pivot position to be close to the position 
of the seed particle [21]. Doing so has been reported to 
extend the operating limit to rf s ~ 0.34. However for 
the case of highly asymmetrical mixtures we find that 
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this strategy does not significantly decrease the number 
of particles in the cluster because as soon as a second big 
particle joins the cluster and is point reflected it causes 
many overlaps with small particles. 



B. Staged insertion algorithm 

Our second MC approach for obtaining effective po- 
tentials for size asymmetrical mixtures is based on the 
staged insertion of a big particle [38-41]. The method 
involves first fixing one big particle at the origin and then 
sampling the free energy change associated with insert- 
ing a second big particle at a prescribed distance rbb from 
the origin. Essentially this amounts to an estimation of 
the chemical potential of the second particle n e x{fbb)- As 
such our approach is close in spirit to one proposed very 
recently by Mladek and Frenkel [42], although their im- 
plementation did not employ staged insertion and was 
therefore restricted to low density systems or those in- 
teracting via very soft potentials. 

The effective potential for big particle separation r^b is 
simply 



W(r bb ) = Hex{rbb) + C 
where the additive constant 



C = - lim Hex{rbb) 

r bb — s-oo 



(2) 



(3) 



can be determined as the excess chemical potential of 
a single big particle in the reservoir of small particles. 
To estimate /J, e x(rbb) we follow the strategy described in 
Ref. [40]. In outline, one imagines that the second big 
particle can exist in one of M possible 'ghost' states in 
which it interacts with a small hard particle (a distance 
rbs away) via the potential 



in) 



m) 



(r6») = -[l-e(2r 6 .-<T6)]]nA< 



Here m = I . . . M (an integer) indexes the ghost states, 
while the associated coupling parameter < A^™) < I 
controls the strength of the repulsion between the big 
particle and the small one. Owing to the step function 
O, the potential is uniformly repulsive over the volume 
of the big particle, and zero elsewhere. Moreover, for 
> o the repulsion is finite so that overlaps between 
small particles and the big one can occur. If we denote by 
N the number of such overlaps at any given time, then 
the configurational energy associated with the ghost big 
particle is 



(4) 



and the big particle is effectively absent from the system. 
To span this range we set the extremal states A^ 1 ) = and 
A< M ) = I, and define some number of intermediate states 
that facilitate efficient MC sampling over the range m = 
1 . . . M, i.e. that permits the ghost particle to fluctuate 
between being a real hard sphere and being absent. The 
measured value of the relative probability of finding the 
system in these extremal states yields the excess chemical 
potential: 



Hex(nb) = In 



p(AW) 
p(A(i)) 



(5) 



Clearly for A^" 1 ) = 0, the big particle acts like a normal 
hard sphere, while for A^ m ' = 1 there is no interaction 



Now since W(rbb) is spherically symmetric, it can be 
estimated from Eqs. 5 and 2, by measuring [i ex (rbb) for 
values of rbb > <Jb along a ID grid. Moreover since each 
such measurement is independent of the others, the ap- 
proach is trivially parallel and thus can be effectively 
farmed out on multi-core processors. 

Details of a suitable Metropolis scheme for sampling 
the full range of m = 1 . . . M have been described in de- 
tail previously [40, 41]. The basic idea is to perform 
grand canonical simulation of the small particles, supple- 
mented by MC updates that allow transitions m — > m± 1 
for the ghost big particle. These transitions are accepted 
or rejected on the basis of the change in the configura- 
tional energy Eq. 4. However, for this strategy to realize 
the aim of sampling the relative probability of the ex- 
tremal states, it is necessary to bias the transitions such 
as to ensure approximately uniform sampling of the M 
ghost states. This is achieved by determining a suitable 
set of weights which appear in the MC acceptance prob- 
ability [43]. 

Additionally it is important to choose sufficient inter- 
mediate states and to place them at appropriate values 
of A such that transitions m — > m ± 1 are approximately 
equally likely in both directions and have a reasonably 
high rate of acceptance. To achieve this we perform a 
preliminary run in which we consider a single big ghost 
particle in the reservoir of small ones. We first define 
M = 1000 values of A in the range (0, 1), evenly spaced 
in In A, and (in short runs) measure the distribution of 
overlaps p(N \X) for each. From this set we then pick 
out those values of A for which successive p(N a \X) exhibit 
an overlap by area of approximately 20%. This criteria 
yields a suitable set of intermediate states. 

Efficiency benefits result from noting that the rate 
of transitions in m depends on how quickly the num- 
ber of overlaps iV relaxes after each successive transi- 
tion. To enhance this relaxation we preferentially per- 
form grand canonical insertions and deletions of small 
particles within a spherical sub- volume of diameter \.2a b 
centred on the second big particle. Updates within the 
sub-volume occur with a frequency 100-fold that outside 
the sub-volume. Our approach -which satisfies detailed 
balance- greatly reduces the time spent updating small 
particles whose coordinates are relatively unimportant 
for the quantity we wish to estimate. 
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A further innovation, applicable to highly size- 
asymmetrical hard sphere mixtures, stems from the ob- 
servation that it is not actually necessary to insert a big 
hard sphere in order to calculate the effective potential. 
Instead it is sufficient and (generally much more efficient) 
to insert a hard shell of infinitesimal thickness. The ba- 
sic idea is that when fully inserted a hard shell parti- 
cle encloses a number of small particles. These remain 
in equilibrium with the reservoir (i.e. their number can 
still fluctuate), but are fully screened from the rest of 
the system because their surfaces cannot penetrate the 
shell wall. Accordingly their contribution to the free en- 
ergy is independent of the position of the shell particle 
and hence their net effect is merely to shift the value 
of the additive constant in the measurement of /U e x (?-&&)■ 
Since the latter is anyway set by hand to ensure that 



lim,. 



W(r 



bb 



0, one does not need to know the 



contribution to the free energy from the enclosed parti- 
cles. 

From a computational standpoint, the task of insert- 
ing a hard shell is much less challenging than that of in- 
serting a hard sphere: essentially the chemical potential 
grows with the particle size ratio like (l/q) 2 rather than 
(1/q) 3 . Consequently, far fewer intermediate stages M 
are required to effect the insertion, which reduces sub- 
stantially the computational expenditure in measuring 
Hex{r) accurately. 

Whilst the staged insertion technique is not as straight- 
forward to implement as the GCA, it does not suffer the 
rapid decrease in efficiency for 77^ > 0.2 observed in highly 
asymmetrical mixtures. It therefore allowed us to attain 
(for q = 0.1), the considerably larger reservoir packing 
fraction of rf a — 0.35. 



C. Correcting for finite size effects 

The effective potential W(r) between two big parti- 
cles is defined in terms of the radial distribution function 
g(r) = gbb(r), with r = r^b, measured in the limit of 
infinite dilution 



-PW(r) = lim ln[ 5 (r)] 

Pb-5-0 



(6) 



for r > <Jb- In our simulation studies this limit is approx- 
imated by placing a single pair of big hard spheres in the 
simulation box. A finite-size estimate to g(r), which we 
shall denote <7z,(r), is then obtained by fixing the first of 
these particles at the origin and measuring (in the form 
of a histogram) the probability of finding the second big 
particle in a shell of radius r — > r + dr, i.e. 



9l{t) 



P(r) 



(7) 



where, the normalization relates to the probability of 
finding an ideal gas particle at this radius: 



Pi 9 (r) = 



Airr 2 

IT 



(8) 



Now the principal source of finite-size error in <7i(r) 
arises from the normalization of Pi g by the system vol- 
ume. Specifically, for a finite sized system, the volume 
occupied by the hard sphere at the origin is inaccessible 
to the second particle. Accordingly, the accessible system 
volume is 



V = V-v 1 



(9) 



where v\ = (1/Q)naf is the hard sphere volume. More 
generally, one should define an effective excluded volume 
V\ for use in Eq. 9, which allows for the fact that the 
small particles can mediate additional repulsions and/or 
attraction between the two big particles. In principle v\ 
is given by 



Vl = 4-7T 



/•oo 

/ [l-ff(r)] 
Jo 



r 2 dr . 



(10) 



It follows from Eqs. 7-10 that the principal finite-size 
contribution to <7i,(r) is just an overall scale factor: 



V 

9(r) = y9L{r) • 



(11) 



Accordingly <7z,(r) approaches V/V at large r instead of 
unity, whilst the calculated effective potential, W(r) de- 
cays to ln(V/V) instead of zero. 

One can conceive of a number of possible approaches 
for dealing with this finite-size error. One is simply to 
minimize it by choosing a very large system volume V so 
that V/V is close to unity. The disadvantage of this ap- 
proach is that in a size asymmetrical mixture, in which 
the big particles are in equilibrium with a reservoir of 
small ones, a huge number of small particles will neces- 
sarily fill the extra space available in a bigger box. All 
the interactions arising from these small particles then 
need to be computed - which can become prohibitively 
expensive. 

Another route, which we have adopted in the present 
work, is to attempt to correct <7z,(r) by estimating the 
overall scale factor in Eq. 11, thus ensuring that g(r) — > 1 
at large r. An expedient approach to doing so, which 
utilizes as much as possible of the information in g L (r) , 
proceeds by determining the cumulative integral of <?i(r): 



G(R) 



f 

Jo 



g L (r)dr . 



(12) 



In practice, this integral was observed to tend towards 
a smooth linear form quite rapidly as the upper limit R 
increases, a fact illustrated for typical data in Fig. 1. A 
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little thought then shows that if when left uncorrected 
g(r) tends to V/V at large r, the limiting gradient m 
of G(R) is m — V/V, which thus provides the requisite 
correction factor for use in Eq. 11. Thus one corrects the 
measured histogram fl'i(r) by first fitting G(R) to obtain 
an estimate of the limiting gradient of the linear part, 
and then scaling gi,(r) according to g(r) = m • 







1 








- g L w 








- G(R) 






Gradient 0.41004 






f 


i , i , i , i 



I I I I I I I I I I I 

1 1.05 1.1 1.15 1.2 1.25 

r/o b , R/o b 



FIG. 1. The measured form of ql{v) for q = 0.05, rf s = 
0.2 obtained for a system size L = 2.5<T(,. Also shown is 
the cumulative integral G(R) = gL{r)dr, together with an 
asymptotic linear fit, the gradient of which yields the finite- 
size correction factor for g(r). 



III. THEORETICAL METHODS 

As mentioned in the Introduction we choose to make 
comparisons between our simulation results and those 
from the DFT insertion method and from the morpho- 
metric and Derjaguin approximations. 



A. DFT 

The DFT insertion method is described in detail in 
Roth et al [14]. It is based on an exact result from 
the potential distribution theorem, for an arbitrary mix- 
ture, that expresses the effective potential between two 
big particles in terms of the one-body direct correlation 
function of the big species in the limit where the den- 
sity pi, of that species vanishes. For DFT treatments of 
hard sphere mixtures that employ the fundamental mea- 
sures theory (FMT) approach [28] the calculation of the 
effective potential requires the computation of the den- 
sity profile p s (r) of the small spheres in the neighbour- 
hood of a single fixed big sphere as well as knowledge 
of the weight functions and the excess free energy den- 
sity of the (binary) mixture [14]. The FMT must be 
sufficiently accurate to describe a very asymmetric bi- 
nary mixture, i.e. one with small values of q, in the limit 



Pb — » 0. In the original DFT studies [14, 18] the Rosen- 
feld (RF) FMT [28] was employed. In the present work 
we employ RF, the White Bear (WB) version [44] and its 
modification the White Bear Mark 2 (WB2) version [45]. 
These versions differ from RF in the choice of the coef- 
ficients 4> a entering the excess free energy function. RF 
yields thermodynamic quantities that are the same as 
those from Percus-Yevick (compressibility) approxima- 
tion, whereas WB incorporates the accurate Mansoori- 
Carnahan-Starling-Leland (MCSL) empirical bulk equa- 
tion of state. In WB2 additional self-consistency require- 
ments are imposed on the pressure. The consistency of 
the WB2 version was demonstrated in calculations of the 
surface tension and other interfacial thermodynamic co- 
efficients for a one-component hard-sphere fluid adsorbed 
at a hard spherical surface [45]. Ref. [46] provides an 
overview of recent developments and describes compar- 
isons between different versions of FMT. 

Botan et al [24] carried out DFT insertion method cal- 
culations as well as explicit (brute force) free energy mini- 
mization for two fixed big spheres using different versions 
of FMT. These authors provide a compendium of the in- 
gredients entering the FMT functionals as well as the 
thermodynamic coefficients required in the morphomet- 
ric approximation and we refer readers to Appendix B 
of [24] for the explicit formulae used in the present cal- 
culations. Their paper is important in pointing to the 
regimes where the DFT insertion method is likely to fail. 
In particular for rf s = 0.419 and q = 0.1 and 0.2 there 
are substantial differences between the results for the de- 
pletion potential calculated by brute force and from the 
insertion method. At higher reservoir packing fractions 
the differences can be even larger. Moreover different 
FMTs can give rise to quite different potentials at high 
small sphere packings. The comparisons made by Oettel 
et al [29] for the depletion force using the RF functional 
suggest that for q = 0.05 the insertion method is not es- 
pecially accurate at rj r s = 0.314 and 0.367. Of course one 
is assuming that the brute force minimization is the more 
accurate method as this requires only a reliable functional 
for a single component hard sphere fluid-not one for the 
asymmetric mixture. 

However, our present study focuses on smaller values of 
rf s than those considered in [24]. Previous studies [14, 18] 
showed generally good agreement between DFT insertion 
results and those of simulation for q = 0.1 and 0.2 and 
rf s typically up to 0.3. Since we are concerned primar- 
ily with investigating the depletion potential for highly 
asymmetrical mixtures in regimes, accessible to simula- 
tion, near the onset of fluid-fluid phase separation, we do 
not concern ourselves with very high values of rf s where 
the DFT insertion method is likely to be inaccurate. 

Another way of viewing this DFT insertion method is 
that it is equivalent [14, 18] to calculating the big-big 
radial distribution function gbb{r) for a binary mixture 
using the test particle route, i.e. one fixes a big sphere 
at the origin and computes the inhomogeneous density 
profile of the big spheres Pb(r) by minimizing the mixture 
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free energy functional for this spherical geometry. Then 
9bb(f;pb) = Pb(r)/pb and the depletion potential is given 

by 



-PW(r)= lim \ng bb {r;p b ) 

Pb^O 



for r > a b . (13) 



In Refs. [14, 18], for all cases considered, it was demon- 
strated that a bulk packing fraction r] b = 10 -4 of the big 
spheres was sufficiently small to ensure that the depletion 
potential calculated from g bb (r) had converged to the lim- 
iting form. In a very recent paper Feng and Chapman 
[47] used the mixture WB theory to calculate g bb (r) via 
the test particle route. For size ratios q — 0.1 they report 
good agreement with existing simulation results [48] for 
concentrations of the big hard spheres as small as 0.002 
and total packing fractions as large as 0.4. However, the 
packing fraction rj b is still too high to be appropriate for 
determining the depletion potential. 

Roth et al [14] also introduced a parameterized form 
for the depletion potential obtained from their DFT in- 
sertion method calculations. Their motivation was to 
provide an explicit form W — 1/2(1/(7 + 1) W(x,Tfl) , 
with x = h/a s and h = r — a b the separation between 
the surfaces of the big spheres, that would be valid for a 
range of size ratios q and reservoir packing fractions 77J 
and therefore efficacious in studies of the phase behaviour 
and correlation functions of binary hard sphere mixtures. 
The authors were influenced by the comprehensive sim- 
ulation studies of Dijkstra et al [4] which had employed 
a very simplified (third order virial expansion ) formula 
for the depletion potential derived in ref. [33]. Roth et 
al aimed to provide a formula, convenient for simula- 
tions of an effective one-component fluid, that captured 
both the short-ranged depletion attraction and the long- 
ranged oscillatory behaviour of W(r). Such a formula 
is of course also useful in making comparisons between 
theory and experimental measurements of the depiction 
potential. Their formula for W consists of a third order 
polynomial at small x and an exponentially damped os- 
cillatory function at large x accounting for the correct 



asymptotic decay [14]. Comparisons made for rf s be- 
tween 0.1 and 0.3 and different values of q showed that 
the parametrized form gave a good fit to the results of 
the numerical calculations. Largo and Wilding [35] em- 
ployed this parametrized form in simulation studies of 
the (metastable) fluid-fluid critical point of the effective 
one-component fluid, comparing their results with those 
from the much simpler parametrized form used in [4]. 

In the present study we noticed that the parametriza- 
tion in [14] did not recover the correct Asakura-Oosawa 
limiting behaviour as nf s — > and this restricts its regime 
of application. Since we are concerned with making di- 
rect quantitative comparisons with the results of our new 
simulations we performed new numerical DFT insertion 
calculations, avoiding parametrization. 
B. Derjaguin and Morphometric Approximations 

A much used theoretical tool of colloid science is the 
Derjaguin approximation [49] that relates the force be- 
tween two large convex bodies immersed in a fluid con- 
sisting of much smaller particles or molecules to the inte- 
gral of the excess pressure of the same fluid contained 
between two parallel walls. In recent years there has 
been considerable discussion about the regime of valid- 
ity of the Derjaguin approximation for our present case 
of a fluid of small hard spheres confined between two 
fixed big hard spheres or between a planar hard wall 
and a single big hard sphere. The reader is referred 
to Refs. [9, 14, 33, 50, 51] should they wish to savour 
the arguments. Herring and Henderson [8, 9] performed 
simulations for the wall-sphere case for q = 0.05 and 
r] r s = 0.3 and 0.4, comparing their results for the deple- 
tion force with those from the Derjaguin approximation 
and from the DFT insertion method [14]. In the present 
work we perform equivalent comparisons, for the sphere- 
sphere case, using what we believe is much more accurate 
simulation data for the depletion potential. 

As shown in [33] the depletion potential difference for 
hard spheres obtained from the Derjaguin approximation 
can be expressed succinctly as: 



£7T , 



W Der (h) - W Der (a s ) = - y((T s + cr b )((7 s - h) 



M)(a s -h) + 2j{ V r a ) 



< h < <7 S 



(14) 



where h is the separation between the surfaces, p{jf s ) is 
the pressure of the small sphere reservoir and j(t)1) is the 
surface tension between a single planar hard wall and the 
small sphere fluid. Within the Derjaguin approximation 
the potential between a wall and a single big sphere is 
precisely twice that between two big spheres: e is 1 for 
sphere-sphere and 2 for wall-sphere. Expressions for the 
pressure and surface tension are listed in Appendix A 
of [24]. Another expression for the surface tension due 



to D. Henderson and Plischke [52] as obtained by fitting 
simulation data was used in [9] . For h > u s the depletion 
potential depends on the excess grand potential of the 
small sphere fluid confined in the planar hard wall slit 
which must be obtained from simulation or DFT [24]. 

Morphometric thermodynamics [30] was developed to 
calculate the solvation free energy (excess grand poten- 
tial) of large convex bodies immersed in a solvent. Its ap- 
plication to determining depletion potentials is described 
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in [24, 29] where it is shown that 



W M or P {h) = -pAV(h)-jAA(h)- KAC(h)-4nR, (15) 

for < h < <r s . Here AV(h) and AA(h) are the volume 
and surface area of the overlap of exclusion (depletion) 
zones around the big spheres (or a wall and a big sphere) 
and AC(h) is the integrated mean curvature of the over- 
lap volume. The thermodynamic coefficients are the pres- 
sure p, surface tension 7 and the two bending rigidities 



k and k ; these four quantities are functions of r/l .The 
fourth term is the difference in integrated Gaussian cur- 
vatures between a dumbbell (47r) and two disconnected 
spheres (8tt) . For h > a s the dumbbell separates into two 
disconnected spheres. Thus WMor P (h > <r s ) = 0. Ex- 
plicit formulae are given in Ref. [24] for the geometrical 
quantities and for the four thermodynamic coefficients. 
Note that WMorpi^j) = —^R{ri r s ), independent of size 
ratio. This term is small in comparison with the others. 

In order to connect with the Derjaguin approximation 
we invoke the colloidal limit , i.e. q — > 0. Then the 
difference 



W M orp(h) - WMorp(Vs) 



-(cr s + a b ){a s - h) 



pW s )(a s -h) + 2 1 {r 1 l) 



K(v r s )K 2 V<<rs ~ h){a s + a b )/2, (16) 



r 



for < h < <r s . Since n is positive the morphometric 
approach contributes an additional attractive term, aug- 
menting the attraction from the pressure (AV(hj) term. 
7 is negative so the surface tension (AA(h)) term gives 
a repulsive contribution to the depletion potential. The 
physical interpretation of the third term in Eq. 15 or 16 
is of a line contribution to the effective interaction asso- 
ciated with the circumference of the edge of the annular 
wedge formed between the two exclusion spheres where 
the line tension is —kit/ 2 [29]. As this term is propor- 
tional to ^fk (not to e ) Derjaguin scaling is violated [24]. 

The morphometric analysis must break down in the 
limit h — > a s where Eq. 15 or 16 predicts that the de- 
pletion force is singular, diverging as (a s — /i) -1 / 2 . The 
reasons for this unphysical limiting behaviour are asso- 
ciated with problems of self-overlapping surfaces as ex- 
plained in Refs. [24, 29]. However, away from this limit 
one might expect the elegant geometrical arguments un- 
derlying the morphometric analysis to capture the essen- 
tial physics. Indeed the comparisons with brute force 
DFT results for the depletion force in Ref. [29] indicated 



rather good agreement for a range of q and rf s = 0.314. 

In Sec. IV we compare the results of Eqs. 14 and 15 
with our simulation data and with results from the DFT 
insertion method. 



IV. RESULTS 

A. Test Case 

We have tested the ability of the GCA to accurately 
determine effective potentials by applying it to the case 
of the Asakura-Oosawa (AO) Model [31, 32]. This 
model describes colloidal hard-spheres in a solvent of 
non-interacting point particles modelling ideal polymer 
that have a hard-particle interaction with the colloids. 
Although not the case of additive hard-spheres which 
is our primary focus in this paper, the extremely non- 
additive AO model does provide a very useful testbed for 
our simulation methodology because the exact form of 
the depletion potential is known, taking the form: [32] 



/3Wao{t) 



0, 



1 - 



or r 

2<r b (l+q) 2aHl+q)3 



<r b < r < a b + a s 



r > a b + a s 



(17) 



where a s is the 'polymer' diameter, i.e. the colloid- 
polymer pair potential is infinite for r < (a b + a s )/2. 

Simulation measurements of g{r) = g bb (r) were per- 
formed for the AO model using the GCA for a system 
comprising a pair of hard spheres in a cubic box of linear 
dimension L = 3a in equilibrium with a reservoir of small 
particles having size ratio q = 0.1. Since the small parti- 
cles are mutually non-interacting, the chemical potential 



of the reservoir is just that of an ideal gas. The depletion 
potential was calculated as j3W(r) = — ln[g(r)] and the 
results were corrected for finite-size effects according to 
the procedure described in Sec. II C. In Fig. 2 we com- 
pare the results of simulations of the effective potential 
with the exact result. Data is shown for various values of 
the reservoir packing fraction. In each instance, the sim- 
ulation results (symbols) are indistinguishable from the 
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FIG. 2. A comparison of GCA simulation measurements of 
the depletion potential of the AO model with the exact ana- 
lytical form. The size ratio is q = 0.1 and data are shown 
for four values of the reservoir packing fraction. Symbols 
are results from GCA measurements of g(r) for a pair of 
big particles, corrected for finite-size effects and transformed 
via f3W(r) = — ln[g(r)]. Statistical errors are smaller than 
the symbol size. Lines are the exact AO effective potential, 
Eq. 17. 



analytical form (lines) within the very small statistical 
errors, a finding that supports the validity and accuracy 
of the simulations. 



B. Effective potentials for additive hard spheres 

We turn now to our measurements of the effective po- 
tential for highly size asymmetrical additive hard spheres 
and the comparison with DFT calculations. Similarly 
to the case of the AO model, our simulations treat the 
small particles grand canonically, i.e. their number fluc- 
tuates under the control of a chemical potential [i r s . The 
value of fjb r s is chosen to yield some nominated value of 
the packing fraction of small particles, rf s , in the notional 
reservoir. Thus the simulations require prior knowledge 
of [i r s {rf s ). In principle, one could employ the Carnahan- 
Starling (CS) approximation [53] to estimate the requi- 
site chemical potential. However in tests we found this 
approximation to be insufficiently accurate for our pur- 
poses. For instance, taking rf s = 0.32 as an example, if we 
employ the CS value for the chemical potential we actu- 
ally measure ff s — 0.3195, which while close to the target, 
lies outside the range of fluctuations in rf s that occur in 
a large simulation box. In order to determine n r s more 
accurately we therefore performed a series of accurate 
grand canonical simulations for the pure fluid of small 
hard spheres in a large box of L = 50a s . We then em- 
ployed histogram reweighting to extrapolate to the pre- 
cise values of the chemical potential that corresponds to 



the various values of rf s that we wished to study. These 
resulting estimates are listed in table I. 



nl 


PA 


0.05 


-1.9079(1) 


0.10 


-0.6770(3) 


0.15 


0.3923(2) 


0.20 


1.5105(1) 


0.32 


5.0472(1) 


0.35 


6.2659(2) 



TABLE I. Measured values of the reduced chemical potential 
Pfis corresponding to each of the packing fractions rf s listed. 
The data were obtained by histogram reweighting the results 
of grand canonical simulations of hard spheres obtained at 
the nearby value of /3/Lts predicted by the CS approximation. 
The simulation cell size was L — 50a s . The definition of fi r a is 
subject to the convention of choosing the thermal wavelength 
to equal the hard sphere diameter. 

Measurements of the radial distribution function g(r) 
were made for a pair of big hard spheres in equilibrium 
with a reservoir of small hard spheres, for the combina- 
tions of values of rf s and size ratio q shown in table II. 
The system size was L = 3<7b for q = 0.1, while for 
q = 0.05, 0.02, 0.01, where the range of the depletion po- 
tential is shorter, L = 2.5<7fc was used. In all cases the 
depletion potential was obtained as j3W(r) = — hx[g(r)] 
with corrections for finite size effects applied as described 
in Sec. IIC. 
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0.05 
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0.01 


0.05 









TABLE II. Combinations of particle size ratio q and reservoir 
volume fraction rf 3 for which we compare simulation estimates 
of depletion potentials with DFT predictions. Values shown in 
normal typeface were studied by simulation using the GCA 
described in Sec. II A, while those in boldface were studied 
using staged insertion MC as described in Sec. II B. 



1. Comparison of simulation and density functional theory 
results 

We now examine a selection of the measured effective 
potentials. Data for q — 0.1, rf s — 0.2 is shown in Fig. 3. 
Despite our use of a rather small histogram bin size of 
just Sr = 0.001 to accumulate estimates of /3W(r), the 
statistical fluctuation is sufficiently small that one can 
simply connect the data points by lines. This allows us to 
better discern differences between the simulation results 
and those of the DFT calculations using the insertion 
method, which are also included on the plot. Data for 
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three versions of DFT are shown, namely the Rosenfeld 
(RF), White Bear (WB) and White Bear 2 (WB2) func- 
tionals. Clearly for these parameters the overall agree- 
ment is very good. To quantify the extent of the accord, 
the two insets to Fig. 3 show a comparison in the range of 
separations close to hard sphere contact (left inset) and 
around the first maximum (right inset). These show that 
near contact, WB2, fares slightly better than WB, which 
is in turn better than RF. Near the first maximum in the 
potential however, the trend is reversed and RF has the 
greatest accord with the simulation data, while WB is 
better than WB2. 




i i i i i i i i i i i 

1 1.1 1.2 1.3 1.4 1.5 

r/o. 

b 



FIG. 3. Simulation and DFT results for the hard sphere de- 
pletion potential /3W(r) for q = 0.1, rf s = 0.2. The abscissa is 
the separation of hard sphere centres expressed in units of the 
big particle diameter Ob- The two insets expand the region 
close to hard sphere contact (left panel) and around the first 
maximum (right panel). 

A similar picture emerges for q = 0.05, rj I = 0.2 as 
shown in Fig. 4. Although here the simulation data is 
not as smooth as for q = 0.1, the form and magnitude of 
the deviations from the DFT are similar. We comment 
later on the results of the morphometric approximation 
Eq. 15. 

Generally speaking, the smaller the size ratio, q, the 
lower the maximum packing fraction rf s for which we can 
obtain good statistics with the GCA. Data for q = 0.02, 
with rf s = 0.1 and rj s = 0.15 is shown in Fig. 5(a) and 
5(b), respectively. For this size ratio and these (low) 
small sphere packings the various versions of FMT per- 
form very well. Data for q = 0.01 with rf s = 0.05 is 
shown in Fig. 6. In this extreme case the insertion DFT 
results are almost indistinguishable from each other and 
from the results of simulation. However one should note 
that there is still a maximum in W(r) - one is not yet in 
the AO limit although the contact value is close to the 
AO value given by Eq. 17. 

For our system, the GCA is operable for rj^ < 0.2. To 
go beyond this limit we have employed the staged inser- 
tion algorithm outlined in Sec. II B. Simulation results 




FIG. 4. As for Fig. 3 but with q = 0.05, n r a = 0.2. Also shown 
are the results of the morphometric approximation Eq. 15. 



for q = 0.1, if s = 0.35 are compared with those from 
DFT calculations in Fig. 7. While the simulation data 
are somewhat noisier, they show that in this regime, quite 
significant discrepancies with the DFT insertion method 
have emerged. The principal form of the discrepancy, i.e. 
DFT underestimates the height of the first maximum, is 
similar in form but greater in magnitude to that seen us- 
ing the GCA at smaller values of rf s (cf. Fig. 3). Once 
again RF fares better than the two WB functionals but 
underestimates the first maximum by about O.bksT. Re- 
sults from the three functionals agree quite well with one 
another and with simulation near contact. 



2. Comparison with Derjaguin and morphometric 
approximations 

In this subsection we make comparisons between our 
simulation and DFT results with those from the approxi- 
mations described in Sec. Ill B. Recall that the Derjaguin 
approximation is specifically designed to tackle small size 
ratios. In Fig. 4 we compare the results of the morpho- 
metric approximation Eq. 15 with those from simulation 
and DFT. Two sets of thermodynamic coefficient were 
used: RF and WB2 [24]. Both versions underestimate 
the maximum of the depletion potential and overestimate 
the magnitude of the potential at contact for q = 0.05 and 
rf 8 = 0.2. By contrast for q — 0.1 and rf s = 0.35, Fig. 7 
shows that both versions of the morphometric approx- 
imation overestimate the maximum and underestimate 
the magnitude of the potential at contact. Fig. 7 also 
shows a pronounced minimum for h close to a s . This 
feature is absent in both simulation and DFT. It is asso- 
ciated with the unphysical divergence of the line tension 
contribution to the depletion force arising in the morpho- 
metric treatment. Recall that WMorph is zero for sepa- 
rations h > er s . 
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FIG. 6. As for Fig. 3 but with q = 0.01, rf a = 0.05. 



FIG. 5. As for Fig. 3 but with (a) q = 0.02, rj r s = 0.1 and (b) 
9 = 0.02,7^=0.15. 



Fig. 8(a) compares the depletion potential difference 
obtained from simulation and insertion method DFT for 
q = 0.1 and 77^ = 0.35 with results from the Derjaguin ap- 
proximation (where plotting the difference is the natural 
choice [33]) and morphometric approximations, Eqs. 14 
and 15 respectively. For this value of q the packing frac- 
tion of the small spheres is sufficiently large to enter the 
regime where fluid-fluid phase separation might occur-as 
discussed below in Sec. IV B 3. Thus it is interesting to 
observe how well these explicit approximations perform. 
Similar remarks apply for q = 0.05 and 77^ = 0.2 for which 
comparisons are presented in Fig. 8(b) 

One sees in Fig. 8(a) that the Derjaguin approximation 
is very poor. Overall the morphometric approximations 
fare considerably better than Derjaguin with RF better 
than WB2 near the maximum. However, both morpho- 
metric versions overestimate the magnitude of the con- 
tact value by about 0.5/csT. Note once again the mini- 
mum close to h/a s — 1 for this packing fraction. The sit- 
uation is clearly different in Fig. 8(b) where the Derjaguin 
and morphometric approximations are reasonably good; 
they bracket the simulation and DFT results. The two 
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FIG. 7. The depletion potential for q = 0.1, rf s = 0.35 
obtained using the staged insertion simulation method and 
DFT. The simulation data points represent the results of 150 
independent measurements of p\V(r) made at various fixed 
values of the big particle separation, though concentrated in 
the range r < 1.12(7;,. Also shown are the results of the mor- 
phometric approximation Eq. 15. 



morphometric versions yield results that are very close 
and even in this difference plot one sees that these fall 
below the simulation results both at maximum and at 
contact. At this smaller value of rf s there is no mini- 
mum visible in the depletion potential. Although plot- 
ting the difference appears to improve the level of agree- 
ment between morphometric and simulation, one should 
recall that it is the actual depletion potential displayed in 
Figs. 4 and 7 which matters, eg. in determining _B 2 (?7s), 
to which we now turn. 
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FIG. 8. Comparison of /3W(h) — f3W{a s ) obtained from sim- 
ulation and DFT for (a) q = 0.1, ^ = 0.35 (cf Fig. 7) and 
(b) q = 0.05, rfs = 0.2 (cf Fig. 4) with results of the Derjaguin 
(Eq. 14) and morphometric approximations (Eq. 15). 



3. Second virial coefficients 

While the various simulation and DFT estimates of ef- 
fective potentials show generally good agreement at low 
rf s , the differences grow with increasing rf s and it is nat- 
ural to enquire as to the likely implications for the prop- 
erties of the bulk mixture and in particular phase be- 
haviour. A useful indicator in this regard is the value of 
the second virial coefficient B 2 : 

B 2 = -2ir (e-^'B^ - l) r 2 dr . (18) 

where the effective pair potential is defined in Eq. 1. 

Previous work by Vliegenthart and Lekkerkerker [54] 
and Noro and Frenkel [55] has shown that an extended 
corresponding states behaviour applies to fluids that 
share the same value of B 2 . Specifically, the measured 
values of B 2 at fluid-fluid criticality were found to be 
similar across a wide range of model potentials. Sub- 
sequent work by Largo and Wilding [35], examined this 



criterion explicitly for the case of two DFT-based hard 
sphere effective potentials which had been fitted to an- 
alytical forms and parameterized in terms of the reser- 
voir packing fraction rf g [4, 14]. Using simulation of a 
single component fluid interacting via a pair potential 
(Eq. 1) with W(r) given by these parameterized deple- 
tion potentials, the value of rf s at which the metastable 
fluid-fluid critical point occurs was determined using an 
accurate approach based on finite-size scaling [35, 56]. 
Interestingly for both q = 0.1 and 0.05 and both choices 
of parameterized potentials, the value of B 2 for the deple- 
tion potential at criticality was in quantitative agreement 
with that of the adhesive hard sphere model (AHS) at its 
fluid-fluid critical point as determined separately in simu- 
lations by Miller and Frenkel [57]. These authors report 
a critical value B 2 HS — -l.207.Bj) where the hard 
spheres second virial coefficient B 2 S = 1-Ku\jZ. The 
level of agreement was much greater than that seen for 
more general model potentials (such as square well or 
Lennard- Jones model studied by Noro and Frenkel) , sug- 
gesting that the quasi-univcrsality of the critical point 
B 2 value holds particularly closely for effective potentials 
whose attractive piece is very short ranged in nature, as 
pertains to highly size asymmetrical hard sphere mix- 
tures. Further confirmation of this has been found very 
recently in simulations of the AO potential where, for 
q = 0.1, Ashton [58] has found that the metastable crit- 
ical point occurs at rf s — 0.249(1), to be compared with 
the prediction 77^ = 0.2482 based on matching to B 2 HS . 

In practical terms the universality of B 2 at the fluid- 
fluid critical point implies that one can predict the critical 
point value of Tf s for effective one-component treatments 
of hard sphere mixtures at small q simply by matching 
the corresponding B 2 to the universal value. Conversely, 
it follows that comparison of B 2 values as a function of rf s 
for different potentials provides a sensitive measure of the 
extent to which their phase behaviour differs. We have 
made such a comparison for effective potentials obtained 
from DFT, the morphometric approximation and simu- 
lation for q = 0.1,0.05 and 0.02. The results are shown 
in Fig. 9(a-c) and demonstrate that at the two larger 
values of q even the relatively small differences that we 
observe between the DFT and simulation estimates of ef- 
fective potentials could lead to significant differences in 
the small particle packing fractions at which fluid-fluid 
phase separation is predicted to occur. Specifically, for 
q = 0.1 based on this B 2 criteria, it seems that the DFT 
with the Rosenfeld functional underestimates the puta- 
tive critical point value of if s by some 13%, while the 
WB2 functional underestimates it by some 9% [59] For 
q = 0.05 (Fig. 9(b)) the discrepancy between DFT and 
simulation has fallen to about 4%, while for q — 0.02 
(Fig. 9(c)), the values of B 2 for the hard spheres mixtures 
arising from the various DFT flavors agree very well with 
one another and with simulation, at least for the range 
of rf s at which bulk phase separation is expected to oc- 
cur. They also agree well with the AO model, suggest- 
ing that the additive and extreme non additive models 
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will have very similar phase behaviour at this value of q. 
Recall that for q < 0.154 the mapping of the binary AO 
mixture to an effective one-component Hamiltonian, with 
the AO depletion potential Eq. 17, is exact and we might 
also expect that for very small q the phase behaviour of 
the full binary hard sphere mixture, at physically rele- 
vant (small) values of rf s , is described accurately by the 
depletion pair potential we calculate here. Many body 
contributions should be negligible. 

Interestingly, the DFT data for B 2 exhibit a broad 
minimum with increasing rf gl as can be seen for q = 0.1 
in Fig. 9(a). The same feature has previously been re- 
ported in Ref . [34] . A similar minimum occurs within the 
DFT for q = 0.05 at rf s m 0.38 (not shown in Fig. 9(b)). 
The origin of the upturn in the value of B 2 beyond the 
minimum appears to be due to the fact that the mag- 
nitude of the first repulsive maximum of W(r) increases 
faster with rf s than the depth of the potential well at 
contact. Unfortunately we could not corroborate the au- 
thenticity of this feature via simulation because it occurs 
at larger values of rf s than are currently accessible to us. 
Should it prove real (rather than being an artifact of the 
DFT), it raises the intriguing possibility that fluid-fluid 
phase separation may occur only within a certain range 
of rf 8 . 

Also plotted in Fig. 9 are the results of the morphome- 
tric approximation Eq. 15 for B 2 . Like the DFT results 
these show minima at all q studied (though only that for 
q = 0.1 is visible in the plotted ranges). For q = 0.1, 
B 2 {rf s ) does not cross the AHS line, implying that the 
theory fails to predict fluid-fluid phase separation at this 
size ratio. At smaller q, the agreement with simulation 
is better, but still poorer than DFT. It is disappointing 
that both versions of the morphometric approximation 
perform poorly for q = 0.02, where we find the results 
are substantially different from those of the AO model. 
The discrepancy with simulation for B 2 appears to arise 
primarily from a failure of the morphometric approxi- 
mation to correctly predict the additive constant in the 
potential, as shown from the comparison of the poten- 
tials of Figs. 4 and 7 with the shifted representation of 
Fig. 8. Whilst morphometric results for the depletion 
force [24, 29] might be in reasonable agreement with DFT 
and simulation, any additive shift is important for B 2 . 



V. DISCUSSION 

In summary, we have employed bespoke MC simula- 
tion techniques to obtain direct and accurate simulation 
measurements of depletion potentials in highly size asym- 
metrical binary mixtures of hard spheres having q < 0.1. 
Small particles were treated grand canonically, the value 
of the chemical potential being chosen to target pre- 
scribed values of the reservoir packing fraction rf s . The 
simulation results were compared with new DFT calcu- 
lations (performed using the insertion method) based on 
the Rosenfeld, White Bear and White Bear Mark 2 func- 
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FIG. 9. (a) Comparison of the second virial coefficient B2 
(Eq. 18) derived from DFT, morphometric and simulation 
measurements of the effective potential for q — 0.1 and var- 
ious rf a . The horizontal dashed line indicates the value of 
b ahs _ _2.527crf for which fluid-fluid criticality was found 
in the adhesive hard sphere model [35, 57]-see text. For val- 
ues of B2 below this line, fluid-fluid phase separation is ex- 
pected, (b) and (c) show corresponding plots for q — 0.05 
and q = 0.02 respectively. Unless error bars are shown, statis- 
tical errors in simulation data points do not exceed the symbol 
size. 



tionals. For < 0.2 generally good agreement with sim- 
ulation was found at all size ratios studied, though on in- 
creasing the packing fraction to 77^ = 0.35 at q = 0.1 sig- 
nificant discrepancies between the various flavors of DFT 
and the simulation estimates were evident. In this latter 
regime, Rosenfeld (RF) was found to be somewhat better 
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than the other functionals in reproducing the height of 
the first maximum of the effective potential, while White 
Bear 2 was marginally the best of the three with regard to 
its prediction for the contact value and for second virial 
coefficients. 

Overall our results show that the DFT insertion 
method provides a reasonably accurate description of 
effective potentials for highly size asymmetrical hard 
sphere mixtures at least in the range of small particle 
packing fractions at which fluid-fluid phase separation 
is likely to occur. Indeed at = 0.35, and q = 0.1 
DFT was found to be more accurate than the morphome- 
tric and Derjaguin approximations, the latter providing 
a particularly poor prediction. This conclusion is partly 
at odds with that of Herring and Henderson [8, 9] who 
assert that both DFT and the Derjaguin approximation 
provide descriptions that are almost equally poor com- 
pared to simulation data (see Fig. 7 of Rcf. [9] which 
refers to q = 0.05 and 77J = 0.3 and 0.4) and advocate 
in particular that neither approach should be used in the 
regime of 'moderate' if s to answer important questions 
such as the existence of fluid-fluid coexistence. While we 
concur with this assessment in the case of the Derjaguin 
approximation, Herring and Henderson's conclusions re- 
garding the accuracy of DFT calculations were reached 
on the basis of simulation estimates for the effective po- 
tential which were not obtained directly, but by integrat- 
ing measurements of the interparticle force as outlined in 
Sec. I A. Perhaps as a consequence, their estimates are 
much noisier (see Fig 6 of [9]) than those presented in the 
present work and consequently -we feel- do not serve as 
a sufficiently reliably indicator of the accuracy of DFT 
especially in the key regime where fluid-fluid phase sep- 
aration might occur. 

Indeed, we have investigated the likely extent of the 
consequences for predictions of phase behaviour arising 
from discrepancies between theory and simulation esti- 
mates of depiction potentials via calculations of the de- 
pendence of the second virial coefficient on rf s . Previous 
simulation studies of phase behaviour in single compo- 
nent fluid interacting via effective potentials [35] have 
shown that when the potential is very short ranged, the 
onset of fluid-fluid phase separation occurs at a near- 
universal value of B 2 = — 2.527of . Based on this crite- 
rion, we found that compared to the effective potentials 
obtained via simulation in the present work, the mor- 
phometric theory provides the poorest predictions of the 
critical packing fraction of small particles (and fails to 
predict phase separation at all at q = 0.1). Those from 
DFT underestimate the critical packing fraction of small 
particles by about 10% for q = 0.1 and about 4% for 
q = 0.05. While these are significant discrepancies, we do 
not feel that they constitute a "qualitative breakdown" of 
the DFT insertion method approach as suggested by Her- 
ring and Henderson [8, 9] on the basis of their simulation 



data, at least not in the regime where phase separation 
is expected. Herring and Henderson speak of a nanocol- 
loidal regime. We interpret this as size ratios q of say 
0.1 to 0.01. Our present study shows that this regime is 
amenable to accurate simulation studies up to values of 
the small sphere packing fraction that are relevant for in- 
vestigations of fluid-fluid phase separation and that DFT 
works well in this regime -the focus of the present paper. 
For larger values of rf s there are serious issues concerning 
the accuracy of the existing DFT approaches and we shed 
no new light on this interesting but somewhat extreme 
regime. 

Turning finally to the outlook for further work on 
highly size asymmetrical mixtures, it would, of course, 
be of great interest to verify the existence (or otherwise) 
of the putative fluid-fluid critical point in the full two 
component size-asymmetric hard sphere mixture. This 
topic remains ebullient. For example, a recent paper 
[60] , based on a version of thermodynamic perturbation 
theory, conjectures that additive hard spheres will ex- 
hibit fluid-fluid separation, albeit metastable with re- 
spect to the fluid-solid transition, for size ratios in the 
range 0.01 < q < 0.1. Our measurements of B2 for the 
depiction potentials obtained in our simulations provide 
(potentially) accurate predictions for the small particle 
packing fraction at which the critical point occurs. We 
arc currently employing a grand canonical version of the 
staged insertion MC method [41] to investigate this mat- 
ter. 

Our simulation methods can be applied to any short 
ranged potential and it would also be of interest to ex- 
amine the influence on the effective potential of adding 
small amounts of finite attraction or repulsion to the bs 
and ss interactions. This would allow us to better model 
real colloidal systems, where one can have a variety of 
interaction potentials, and where, even in systems (such 
as sterically stabilized PMMA) which approximate hard 
spheres rather well, one expects residual non-hard sphere 
interactions [61]. Previous work [34, 62] has suggested 
that the effects of such residual interactions may be rep- 
resented in terms of a non-additive hard sphere mixture. 
It would be of interest to examine this proposal explicitly 
using accurate simulation data and DFT calculations. 



ACKNOWLEDGMENTS 

This work was supported by EPSRC grant 
EP/F047800 (to NBW). Simulation results were 
partly produced on a machine funded by HEFCE's 
Strategic Research Infrastructure fund. The authors are 
grateful to Martin Oettel for comments on the original 
manuscript and in particular on the morphometric 
approximation. 



16 



C. N. Likos, Phys. Rep. 348, 267 (2001). 
W. G. McMillan and J. E. Mayer, J. Chem. Phys. 13, 
276 (1945). 

L. Belloni, J. Phys. Condens. Matter 12, R549 (2000). 
M. Dijkstra, R. van Roij, and R. Evans, Phys. Rev. E 
59, 5744 (1999). 

P. Germain and S. Amokrane, Phys. Rev. Lett. 102, 
058301 (2009). 

W. C. K. Poon, A. D. Pirie, and P. N. Pusey, Faraday 
Discuss. 101, 65 (1995). 

P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, 
F. Sciortino, and D. A. Weitz, Nature 453, 499 (2008). 

A. R. Herring and J. R. Henderson, Phys. Rev. Lett. 97, 
148302 (2006). 

[9] A. R. Herring and J. R. Henderson, Phys. Rev. E 75, 
011402 (2007). 

H. N. W. Lekkcrkcrker and R. Tuinier, Colloids and the 
Depletion Interactions, Lecture Notes in Physics, Vol. 833 
(Springer Berlin / Heidelberg, 2011). 

D. Rudhardt, C. Bechinger, and P. Leiderer, Phys. Rev. 
Lett. 81, 1330 (1998). 

C. Bechinger, D. Rudhardt, P. Leiderer, R. Roth, and 
S. Dietrich, Phys. Rev. Lett. 83, 3960 (1999). 
J. C. Crocker, J. A. Matteo, A. D. Dinsmore, and A. G. 
Yodh, Phys. Rev. Lett. 82, 4352 (1999). 
R. Roth, R. Evans, and S. Dietrich, Phys. Rev. E 62, 
5360 (2000). 

D. Kleshchanok, R. Tuinier, and P. R. Lang, J. Phys: 
Condens. Matter 20, 073101 (2008). 
T. Biben, P. Bladon, and D. Frenkel, J. Phys: Condens. 
Matter 8 (1996). 

R. Dickman, P. Attard, and V. Simonian, J. Chem. Phys. 
107, 205 (1997). 

B. Gotzelmann, R. Roth, S. Dietrich, M. Dijkstra, and 
R. Evans, Europhys. Lett. 47, 398 (1999). 
P. Attard, J. Chem. Phys. 91, 3083 (1989). 
J. Malherbe and S. Amokrane, Mol. Phys. 99, 355 (2001). 
J. Liu and E. Luijten, Phys. Rev. E 71, 066701 (2005). 
S. A. Barr and E. Luijten, Langmuir 22, 7152 (2006). 
Though see [18] for a grand canonical study at q — 0.2. 
V. Bo^an, F. Pesth, T. Schilling, and M. Oettel, Phys. 
Rev. E 79, 061402 (2009). 

S. B. Yuste, A. Santos, and M. L. de Haro, J. Chem. 
Phys. 128, 134507 (2008). 

D. Goulding and S. Melchionna, Phys. Rev. E 64, 011403 
(2001). 

D. Goulding, Ph.D. thesis, University of Cambridge 
(2000). 

Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989). 
M. Oettel, H. Hansen-Goos, P. Bryk, and R. Roth, Euro. 
Phys. Lett. 85, 36003 (2009). 

P.-M. K6nig, R. Roth, and K. R. Mecke, Phys. Rev. 
Lett. 93, 160601 (2004). 

S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 
(1954). 



[32] S. Asakura and F. Oosawa, J. Polym. Sci. 33, 183 (1958). 
[33] B. Gotzelmann, R. Evans, and S. Dietrich, Phys. Rev. 

E 57, 6785 (1998). 
[34] R. Roth, R. Evans, and A. A. Louis, Phys. Rev. E 64, 

051202 (2001). 

[35] J. Largo and N. B. Wilding, Phys. Rev. E 73, 036115 
(2006). 

[36] C. Dress and W. Krauth, J. Phys. A 28 (1995). 
[37] J. Liu and E. Luijten, Phys. Rev. Lett. 92, 035504 (2004). 
[38] I. Nezbeda and J. Kolafa, Mol. Sim. 5, 391 (1991). 
[39] P. Attard, J. Chem. Phys. 98, 2225 (1993). 
[40] N. B. Wilding and M. Muller, J. Chem. Phys. 101, 4324 
(1994). 

[41] D. J. Ashton and N. B. Wilding, Mol. Phys. 109, 999 
(2011). 

[42] B. M. Mladek and D. Frenkel, Soft Matter 7, 1450 (2011). 
[43] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, 

and P. N. Vorontsov-Velyaminov, J. Chem. Phys. 96, 

1776 (1992). 

[44] R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: 

Condens. Matter 14, 12063 (2002). 
[45] H. Hansen-Goos and R. Roth, J. Phys. Condens. Matter 

18, 8413 (2006). 
[46] R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010). 
[47] Z. Feng and W. G. Chapman, Mol. Phys. 109 (2011). 
[48] M. Alawneh and D. Henderson, Mol. Phys. 106, 607 

(2008). 

[49] B. V. Derjaguin, Kolloid-Z 69, 155 (1934). 
[50] J. R. Henderson, Physica A 313, 321 (2002). 
[51] M. Oettel, Phys. Rev. E 69, 041404 (2004). 
[52] D. Henderson and M. Plischke, Proc. Roy. Soc. A 410, 
409 (1987). 

[53] J. P. Hansen and I. R. McDonald, Theory of Simple Liq- 
uids (Academic, London, 2006). 

[54] G. A. Vliegenthart and H. N. W. Lekkerkerker, J. Chem. 
Phys. 112, 5364 (2000). 

[55] M. G. Noro and D. Frenkel, J. Chem. Phys. 113, 2941 
(2000). 

[56] N. B. Wilding, Phys. Rev. E 52, 602 (1995). 
[57] M. A. Miller and D. Frenkel, J. Chem. Phys. 121, 535 
(2004). 

[58] D. J. Ashton, Unpublished results. 

[59] We note that in Refs. [34, 54] an empirical (average) value 
gcrit _ _j 5gHS wag usec j £ estimate the critical point. 
We prefer the AHS value as an indicator of the onset of 
phase separation since we focus on short-ranged (sticky) 
potentials, following Largo and Wilding [35]. 

[60] P. Sillren and J.-P. Hansen, Mol. Phys. 108, 97 (2010). 

[61] P. Germain, J. Chem. Phys. 133, 044905 (2010). 

[62] A. A. Louis and R. Roth, J. Phys.: Condens. Matter 13, 
L777 (2001). 



